A statistical mechanics approach to the factorization problem 



O 

(N 

D 
tin 



O 

CD 



O 
o 



On 
O 



X 



p. Henelius^ and S. Girvin'^ 

^Theoretical Physics, Royal Institute of Technology, SE-106 91 Stockholm, Sweden 
^Department of Physics, Yale University, P. O. Box 208120, New Haven, CT 06520-8120 

(Dated: January 11, 2013) 

We map the problem of finding the prime factors of an integer to the statistical physics problem 
of finding the ground state of a long-range Ising-like model. As in the strongly disordered Newman- 
Stein (NS) spin-glass model, the bond distribution is exponentially wide and grows with system size, 
but unlike the NS model we find that it is not wide enough for a greedy algorithm to be applicable. 
On the other hand, we also find that the frustration and exponential width of the bond distribution 
renders classical and quantum annealing and tempering methods no faster than a random search 
for this challenging model. 

PACS numbers: 75.10.Hk,75.10.Nr,75.10. Jm,75.40.Mg 



I. INTRODUCTION 

Since the security of public-key cryptography de- 
pends on the presumed difficulty of computing the fac- 
tors of large integers, there is very strong incentive to 
find efficient factorization algorithms. Shor's quantum 
algorithnJi which can find the factors in a time that scales 
polynomially with the system size has motivated a mas- 
sive research effort into quantum computing. While there 
are no known classical algorithms that scale polynomi- 
ally, the state-of-the-art general number field sieve scales 
sub-exponentiallyi^ The time required to factor a com- 
posite integer q scales as exp[(log(7)3 (loglogg)^], and 
this has enabled factoring of, for example, the RSA-200 
(200 digits, 663 bits) challenge in 2005. 

Our aim in this paper is to view the problem from a sta- 
tistical physics point of view. We show that the problem 
of finding two prime factors of a large composite integer 
can be transformed to an optimization or statistical me- 
chanics problem in which Ising spins are used to represent 
the numbers in binary format. The resulting expression 
can be treated as a system of interacting two-state spins, 
similar to an Ising model in statistical physics with frus- 
trated long-range multi-spin interactions. The model can 
be formulated so that the desired prime factors are given 
by the lowest energy configurationi^i^ 

For statistical mechanics problems without frustration 
and with moderate or no disorder, a number of very effi- 
cient Monte Carlo algorithms have been developed. On 
the other hand, certain special classes of disordered opti- 
mization problems, such as the minimum spanning tree 
problem,— are exactly solved by simple "greedy" algo- 
rithms, where locally optimal choices lead to a global so- 
lution. Somewhat paradoxically, a broader class of prob- 
lems is approximately soluble by greedy algorithms in 
the limit of extreme disorder. For example, the resis- 
tance of a random resistor network^ can be solved to a 
good approximation by selecting the single cheapest path 
that percolates across the cluster. This is asymptoti- 
cally exact in the limit of infinite disorder. Finding the 
ground state of certain models of classical spin glasses^ 



are also minimum spanning tree problems in the limit of 
very broad disorder. For quantum spin problems, asymp- 
totically exact renormalization group methods based on 
analogous ideas have been developediSdi. Here the pair 
of spins with the single strongest bond is integrated out 
which perturbatively modifies the remaining bonds. As 
long as the bond distribution continues to broaden under 
repetition of this renormalization procedure, the system 
flows to an infinite-disorder fixed point and this treat- 
ment is asymptotically exact. 

As we will see below, the factoring problem presents 
several extreme difficulties. First we require the exact 
ground state. Nearby solutions are of no use. Sec- 
ond, the problem contains multi-spin interactions and the 
bond strengths, while not truly random as in a spin-glass 
model, are frustrated and very broadly distributed, but 
not so broadly that a greedy algorithm is even approx- 
imately applicable. Third, the energy landscape is ex- 
tremely rough and there is no concept of nearness in spin 
space. Two nearly identical spin configurations can have 
vastly different energies and conversely two very different 
configurations can have nearly identical energies. 

Our main result is the formulation of the problem in 
a spin language which sheds light on why this is such 
a difficult optimization problem. We explore a number 
of numerical methods for attacking this model, including 
greedy algorithms, single and multi-spin updates, paral- 
lel tempering and quantum annealing. However we find 
that all exhibit exponential cost and, being essentially no 
faster than a random search, are completely defeated by 
this difficult model. 



II. THE ISING MODEL 

In this section we discuss the form of the two-state 
model representing the factorization problem. Assume 
that we want to find two (unknown) prime factors pi 
and p2 of a composite integer q. We can write the factors 
in binary form, pi = YJ^Iq T^S} and p2 = YTjZl 2-'S'J, 
where the spin variables S € {0, 1}. Considering two n- 
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spin factors, the factorization problem, q 
written 



n-l 
hj=0 



P1P2, can be 



(1) 



The right hand side is of the same form as a long-range 
Ising model, H = j JijS} , with the exponential 
coupling constant Jij — 2*+^ . Note that only spins be- 
longing to separate integers interact. We can now use 
Eq. ([T]) to construct a cost function that can be used as 
an effective Hamiltonian in a Monte Carlo simulation. 
The ground state should be given by the desired prime 
factors and first we consider 



n 


Pi 


P2 


q 


10 


601 


911 


547511 


12 


2081 


3329 


6927649 


14 


10007 


15091 


151015637 


16 


40093 


60013 


2406101209 


18 


150011 


140007 


36003690077 


20 


700057 


900001 


630052000057 


22 


2500339 


3500227 


8751754076953 


24 


11600489 


14000083 


162407808840587 


26 


41615281 


61616479 


2564187087815599 


28 


150243361 


220293523 


33090127134000803 


30 


800000087 


900000083 


720000144700007221 



TABLE I: Pregenerated composite integers along with factors. 



III. GREEDY ALGORITHMS 



H„ 



(2) 



where m is a positive number. The ground state is doubly 
degenerate due to a simple permutation of pi and p2 ■ We 
will consider the cases m = 1/2,1,2. The case m = 2 
has the simplest interpretation as a statistical mechanics 
problem since it corresponds to a model with two- and 
four-spin interactions. A similar model has been studied 
recently in the context of a quantum adiabatic algorithm^ 
and a branch and bond optimization methodi^ 

Notice for this problem that the bond strengths are 
powers of two and therefore cover an enormous range. 
On the other hand, we also note that the strongest 
bond approximately equals the sum of the other bonds, 
127=0 2' = 2" — 1, a limit to which we will return in the 
section on greedy algorithms. Flipping all other spins in 
an integer can therefore approximately cancel the effect 
of flipping the highest spin. As a simple example consider 
the drastic effect of " carrying" in the following addition 
01111111 -Kl = 10000000, which clearly demonstrates the 
lack of a simple connection between Hamming distance 
in spin space and the distance in energy. 

Introducing a fictitious temperature T — 1/ f5 the ther- 
mal expectation value of the energy is given by (H) = 
J2i EiS''^^' /J2i e~^"^S where the sum is over all states, 
and Ei is the energy of state i. It is worth noting that the 
factorization problem differs from the spin glass problem 
in that we know the lowest energy eigenvalue, namely 
zero, and we need to find only the lowest energy eigen- 
state. During a stochastic simulation one can therefore 
immediately interrupt the execution of the program when 
the true ground state is found, which makes the factor- 
ization problem an ideal testing case for ground state 
algorithms. We have generated several instances of the 
factorization problem for increasing system size, and in 
TableUwe list the series of composite integers, along with 
the prime factors, used in this work. 



The Edward-Anderson model, H = JijSiSj, is 

the archetypical Ising spin glass model. The coupling 
parameters Jy are quenched random variables, for exam- 
ple Gaussian distributed with mean zero. The combina- 
tion of disorder and frustration makes spin glass models 
very challenging to solve, and there is no known gen- 
eral solution. However, in the limit of very broadly 
distributed coupling constants Jij finding the ground 
state of this model maps to a minimum spanning tree 
problem, solvable by a greedy algorithmj^i^ The crite- 
rion is that the magnitude of each coupling is greater 
than the absolute sum of all smaller couplings, which, 
as a consequence, means that the width of the distri- 
bution increases with system size. This may, at first, 
seem to hold also for the factorization model, Eq. ([1]), 
since J27=o 2* = 2" — 1 < 2". However, there are two 
caveats. First, all spin pairs for which i+j = k share the 
same coupling constant Jk, so there are many coupling 
constants with the same magnitude. Second, when the 
model is formulated so that the prime factors are given 
by the ground state, like H2, the resulting model contains 
multi-spin interactions. 

Next we consider the limit of an even broader bond 
distribution and show that it leads to important simpli- 
fications also for the factorization model. As mentioned 
above there are, in general, many couplings of a given 
magnitude Jk- They all contribute to the fc:th digit of 
the composite integer q. In addition a carry bit must be 
added. If we can avoid the carry bit the problem simpli- 
fies in that the bonds of magnitude Jk directly determine 
the fc:th digit of q. This means that we can, in principle, 
satisfy the bonds (and corresponding digits in q) in de- 
creasing order without running the risk of weaker bonds 
upsetting already satisfied stronger bonds. 

As an illustration, we consider prime factors of the 
form p = 'E'l^ok'S^, where S G {0,1}, but k > 2. If 
k > n no carry operation is needed, and here we consider 
the case of fc = 10. Examples of such prime factors in- 
clude (in base 10 representation): 11, 101, 1011, 101111, 
1011001, 1100101. The simplest case is 11 x 101 = 1111. 
Using long multiplication we work out the factors simul- 
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taneously from the high (since there are no carry bits) 
and low end: 

c b a 
X fed 
cd bd ad 
ce be ae 
cf bf af 
1111 

At the high end cf — and ce + bf = 1 imply that c=l, 
f=0 (or vice versa) and e=l. At the low end ad = I im- 
plies that a = d = 1. Both remaining conditions yield 
b=l rendering the factors 11 and 101. Solving the prob- 
lem decade by decade is therefore much simpler in the 
limit of a bond distribution broad enough to prevent the 
carry operation. Thus, while the distribution of bond 
strengths in the full factorization problem is exponen- 
tially broad, it is not broad enough to allow us to apply 
a greedy algorithm in which we satisfy the bonds one at 
a time. 



IV. STOCHASTIC ALGORITHMS 

Stochastic Monte Carlo methods have long been used 
to calculate thermal expectation values, as well as ground 
state properties, of statistical models. Depending on the 
model, and method, large systems can often be studied 
to high precision. Using cluster algorithms the critical 
temperature and other properties for three-dimensional 
Ising-like models have been determined up to 6 deci- 
mals by performing simulations of systems containing 
128'^ spinsi^, and the ground state energy of the two- 
dimensional Heisenberg model has been determined to 
5 decimals using quantum Monte Carlo methods^ii^. 
There are also many examples in which the use of ef- 
ficient algorithms can change the functional dependence 
of the scaling with respect to computational effort. The 
use of cluster updates at the critical point for the Ising 
model in effect eliminates the problem of critical slowing 
down, and the computational effort decreases from 
for single-spin updates to L'^ for cluster updates, where d 
is the dimensionality of the model and L is the linear sys- 
tem size In some cases use of the loop algorithms can 
make the exponentially difficult sign problem tractable 
in polynomial timeJ^^ 

This suggests that Monte Carlo methods developed for 
statistical physics problems could be useful in analyzing 
the factorization problem, since finding factors contain- 
ing a few hundred bits is considered a hard problemi^ 
However, the Ising model that arises from the factoriza- 
tion problem has a very complex energy landscape, with 
multiple local minima separated by very high barriers. 
Factorization is but one of many such difficult optimiza- 
tion problems whose solution requires finding a global 
minimum in a very complex landscape. Other examples 
include various aspects of circuit design in electronics,^^ 



protein folding in life science^, spin-glass behavior in ma- 
terials science2£ and the traveling salesman problerri^ii^ 
in computer science and mathematics. In an attempt to 
alleviate the problems associated with the complicated 
energy landscape the methods of thermal annealing^"^ and 
parallel tempering^* have been developed. A fictitious or 
real temperature is introduced and as this parameter is 
lowered, the system settles in a local minimum. If the 
cooling is sufficiently slow the ground state is found<^ 
However, for many complex systems it is practically im- 
possible to reduce the temperature slowly enough. 

Classical simulated annealing relies on thermal fluctu- 
ations to find the ground state. If the energy barriers 
separating different minima are high and narrow, quan- 
tum tunneling can be a more efficient way to equilibrate. 
Quantum mechanical systems are able to tunnel through 
barriers, instead of traversing the barriers. In quan- 
tum annealing the minimum is found using the quan- 
tum mechanical tunneling effect instead of thermal fluc- 
tuations. The efficiency of quantum annealing has been 
studied both experimentally^^ and computationally^Zr— , 
in which case quantum annealing can be realized by intro- 
ducing off-diagonal terms into the classical model, which 
cause tunneling between the diagonal, classical states. 

Next we consider a number of different approaches to 
a Monte Carlo simulation of the factorization problem in 
order of increasing complexity. 



A. Random search 

The simplest stochastic method for solving the factor- 
ization problem is to generate random integers p and 
interrupt the search when qmodp = 0. The probabil- 
ity of finding a factor is fj, = 2/2" for a single attempt 
(since there are two factors). Since the probability is con- 
stant for each attempt the number of trials before success, 
P{N), is Poisson distributed, P{N) — /iexp(— /lA^), with 
an average of iV = /i^^ = 2"^^. 

We use the random search as a point of reference and 
next we consider Monte Carlo techniques using impor- 
tance sampling. If the weight function for a problem 
is fairly smooth and dominated by a few pronounced 
maxima, then it may be possible to generate states dis- 
tributed according to the weight function, with a sub- 
stantial gain in performance. However, the weight func- 
tion for the factorization problem is very complex and the 
potential gain less certain. To demonstrate the compli- 
cated structure of the factorization problem we plot the 
function q modp for the case of g = 547511 — 601 x 911 
in Fig. [1] We note an average linear increase of q mod p, 
but otherwise there is no apparent ordered structure. In 
the next sections we investigate whether importance sam- 
pling can nevertheless still be used to improve conver- 
gence to the ground state. 
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FIG. 1: The remainder after dividing q=547511 by p. The 
remainder is equal to zero for p= 601 and 911 



B. Simple spin flips and local temperature 
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FIG. 2: The average number of attempts, N before finding a 
prime factor as a function of temperature. The temperature 
is shown in units of 2^". The curves are for n=22,18,14 and 
10-spin factors from top to bottom. 



Using the Boltzmann weight for state i, Wi = 
exp{—l3Ei), we implement a single-spin flip Metropolis 
algorithm for the model defined by Eq. ([2]). As a a sin- 
gle spin flip is attempted the new weight function W' 
is calculated, and the spin flip is accepted with proba- 
bility p = ma.x{W' /W, 1) = max(exp(— /3Ai?), 1), where 
AE — E' — E. To investigate the performance of the 
Metropolis algorithm we start with two randomly chosen 
trial factors pi and p2 ■ During the simulation we choose 
spins at random and attempt to flip single spins with the 
above probability. When one of the factors is found the 
execution is halted, and the number of attempts, N is 
recorded. Note that all the states visited during the sim- 
ulation are counted, independently of whether the state is 
accepted or not. This is repeated until a reliable average, 
N can be calculated (typically one thousand runs). We 
restrict the search to odd factors, by locking the lowest 
spin in the 1-state. This is to prevent one of the factors 
becoming zero, in which case the energy is independent 
of the value of the other factor. In Fig. [2] we display the 
results obtained for some of the integers listed in Tab. HI 

For each system size there is a minimum at an inter- 
mediate temperature, and as the temperature is lowered 
the number of attempts increases dramatically, while it 
levels out at as the temperature is increased. This be- 
havior is understood considering Fig.[3l where the accep- 
tance probability for the different spins are displayed for 
the five temperatures recorded in Fig. [5] for n=22. At 
the highest temperature the acceptance rate approaches 
unity for all the spins. This means that the spins are 
freely fiuctuating, and the results approach the random 
search described above (consecutive configurations are 
still correlated, so the random search is faster than the 
high temperature limit of the Metropolis algorithm) . As 
the temperature is lowered, the high spins feel the effect 
first. For the nearest neighbor Ising model the change 
in energy, AE = E' — E is oi the order J, the uniform 




spin index 



FIG. 3: The acceptance rates for attempted spin flips of the 
22-spin 



coupling constant, for all spins considered. However, for 
the model described by Eq. ([2]), the change in energy is 
of order 2^ x p2 when an attempt is made to flip the j:th 
spin in pi. As the temperature is lowered, the proba- 
bility of flipping the high spins decreases quickly. The 
minimum in Fig. [5] occurs when the probability of flip- 
ping the highest spin is about 0.05, just before it freezes. 
Lowering the temperature further causes the number of 
attempts to increase dramatically. 

This indicates that one could adjust the model, given 
by Eq. ([2]) so that the probability distribution for accept- 
ing a spin flip is more even. Any alterations are allowed 
as long as the ground state is unchanged, and the aim is 
to decrease the level spacing at higher energies. There- 
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fore we consider the following forms of the Hamiltonian, 



Hi 



'2'^' sis': 



(3) 



and 



Fin = In 



-J2T+^SlS^ 



(4) 



Both the square root and the logarithm are monotoni- 
cally increasing functions that do not change the order 
of the states. The logarithm function, in particular, ap- 
proximately cancels the exponential factor 2^ x p2 and 
allows all the spin to be updated in a more even fashion. 

Yet another way to increase the fluctuations of the 
higher spins is to introduce a site-dependent tempera- 
ture, Ti. A higher temperature at the higher spins en- 
sures more even fluctuations. We therefore define a local 
temperature Ti ^ T xk"^, with a parameter k whose value 
can vary from 1 (no change) to 2 (cancels the exponen- 
tial factor 2^ x ^2)- In Fig- HI we compare the efficiency 
of the different approaches. We display the number of 
states visited before a prime factor is found as a func- 
tion of system size. The temperature is adjusted for each 
data point to optimize the search. All methods are com- 
pared to the random search, for which N = 2"~^, since 
we restrict the search to odd integers, and there are two 
distinct prime factors. We notice that all methods rep- 
resent slight improvements over the random search, but 
no model works better than Hi (absolute value), which 
scales like 0.58 x 2""^ and requires about half the number 
of visited states compared to a random search. 
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FIG. 4: The number of states visited before the n-spin prime 
factors are found for a random search. Hi (absolute value), 
Hi (square root), Hin (logarithm) and local temperature Ti = 

T X 1.5\ 

In Fig. [S] we compare the acceptance rates for the dif- 
ferent methods. The acceptance rate for Hi (absolute 



value) and Hi (square root) quickly saturate to unity, 
while H\n (logarithm) and a local temperature scaling 
like Tj = T X 1.5' result in a much more even accep- 
tance rate. The results clearly show that obtaining a 
more even acceptance rate does not necessarily speed up 
convergence to the ground state. 




10 

spin index 



FIG. 5: The acceptance rate for spin flips as a function of 
spin index for a 22-spin factor. Rates displayed for Hi (ab- 
solute value). Hi (square root), Hin (logarithm) and local 

temperature Ti = T x 1.5\ 



C. Parallel tempering 

In the previous section the behavior of the model was 
investigated while the temperature was held constant. 
More efficient algorithms for converging to the ground 
state rely on a temperature that evolves during the sim- 
ulation. The main methods are annealing22jSS,, where the 
temperature is decreased as the simulation progresses, 
and tempering methodsSi, where the temperature fluctu- 
ates between a maximum and a minimum during the sim- 
ulation. Next we implement a parallel tempering method 
for the factorization problem. In the factorization prob- 
lem the energy scale and the position of the spins are 
closely tied together. Once the energy has been lowered 
sufficiently, the higher spins are entirely frozen, and if 
they are not in the correct positions the ground state 
cannot be found. Compared to annealing methods the 
tempering method offers the advantage that the system 
can return to higher temperatures and explore several 
local minima. 

In the method of parallel tempering several copies, or 
replicas, of the system are simulated concurrently. Each 



replica is initially assigned a temperature, Ti, and after 
performing a number of Monte Carlo updates at the as- 
signed temperatures attempts are made to swap nearby 
temperatures with a probability P(Ti,Tj) = eyip{Ei — 
Ej)/{Ti — Tj), which preserves detailed balance. The 
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attempted temperature swaps are repeated at regular in- 
tervals, and in this manner the temperature of a given 
replica varies during the simulation. The method has 
been very successful in equilibrating disordered spin-glass 
systems at low temperatures as well as studying phase 
transitions. Given the high energy barriers between low 
lying states of the factorization model one could expect 
parallel tempering to allow the replicas to transverse the 
barriers and not get stuck in a given local minimum so 
easily. 

We have implemented a parallel tempering algorithm 
for the factorization problem, based on the Metropolis al- 
gorithm of Hi. One attempt is made, on average, to flip 
every spin in all the replicas, and thereafter an attempt is 
made to switch all neighboring temperatures. The max- 
imum temperature is set to where the Metropolis accep- 
tance rate for flipping the highest spin is about 0.9, and 
the lowest temperature when the acceptance rate for the 
lowest spin is about 0.1. In Fig. IH] we show the tem- 
perature fluctuations for a single replica of a n=40-spin 
factorization problem with pi = 1000000000003, p2 = 
500000000023, g = pi * = 50000000003800000006069. 
The temperature varies from a maximum of T^ax — 
1 X lO^'' to a minimum of Tmin = 1 x 10'^"'^. In between 
there are 40 temperatures that ensure that the swap rates 
remain above 0.5. As can be seen from the figure the 
temperature of the replica wanders repeatedly back and 
forth between the highest and lowest temperatures. 




2000 4000 6000 8000 
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FIG. 6: The temperature fluctuations for replica during par- 
allel tempering 



However, counting the total number of attempted spin 
updates in all the replicas until a prime factor is encoun- 
tered we find that the method is not more efficient than 
the Metropolis algorithm considered in the last section. 
This indicates that although the tempering method im- 
proves equilibration at low and finite temperatures, it 
is not, in this case, superior in picking out the ground 
state itself. There are many low-lying energy states, and 
finding precisely the right one is a difficult problem. 



D. Classical SSE cluster update 

So far we have discussed single spin flips, but one goal 
of this investigation is to determine whether cluster up- 
dates, which have proved highly useful for many difficult 
problems, can be of use for the factorization problem. 
Since flipping the j:th spin in factor pi changes the en- 
ergy by 2^ X p2 one could, in principle, offset the large 
energy change by simultaneously flipping several spins. 
Cluster updates for long-range classical models have been 
developed within the framework of the Swendsen-Wang 
update^^ and the stochastic series expansion'^'*. How- 
ever, the factorization models lacks the up-down symme- 
try of the standard Ising model, which cluster updates 
usually rely on. If the spin variables assume the val- 
ues ±1 the product SiSj is unchanged if both spins are 
flipped. This is not the case for the factorization model 
of Eq. since S takes the values and 1. By a transfor- 
mation S' = 2S — 1 we can introduce a new variable S' 
that takes the values ±1, but this introduces single spin 
operators S'^ in Eq. which also destroy the up-down 
symmetry. Therefore we have not been able to introduce 
a large-scale cluster update, but instead we implement 
a "small-cluster" update within the SSE method, which 
we describe next. 

The SSE method is based on a Taylor expansion of the 
partition function Z, 

a n—Q 

where \a) is a complete set of basis states. The SSE 
method has been applied to Ising model with arbitrary 
interactions, and we refer to Ref. i34j for a detailed de- 
scription. Here we only describe modifications that arise 
when applying the method to the factorization model. 

Expressing the multiplication of two integers in binary 
form we obtain the model 

n 

^ = E J^^J^lS], (5) 

with Ji j — 2'+-'. Dcflning the bond operator Hij = 
Jij{l ~ SlSf) this can be written as 

n n 

^= E -^■^+ E (6) 

Including additional unit operators /, the Taylor expan- 
sion can be written as 

a n=0 Sl 

where we have introduced a cut-off at n = i, and Sl is 
the operator string Sl — Ylp=iHp with Hp £ {Hij,I}. 



The matrix element in Eq. ([7]) can be written as a product 
of elements of the form 

{SlSj\Hij\SlS'^), 

and for the model we consider only the matrix elements 
{00\H,j\00) = (01|if^j|01) = (10|i?,j|10) = J^j con- 
tribute since {ll\Hij\ll) — 0. 

In order to sample the configuration space of all opera- 
tor sequences Sl and all states |a), two types of updates 
are necessary. The first update changes the number of 
non-identity operators in the sequence by attempting to 
exchange identity and bond operators. The probability 
of exchanging a unit operator for a bond operator is 



Phr„ 



bond 



and the reverse probability of exchanging a bond opera- 
tor for a unit operator is 



Punit — 



L 



1 



If a bond operator is to be inserted, the particular bond 
is chosen according to the relative weight of the bond, 
as described in Ref. HI- 

The update described above only changes the opera- 
tor sequence Sl, and the state \a) is not affected. The 
classical cluster update described in Ref. [s^ flips all the 
spins that are interconnected by bond operators; since 
the Ising operators depend on the relative orientation of 
the spins this does not change the weight and is always 
allowed. For the factorization model flipping two spins 
in the state |0) connected by a bond leads to the forbid- 
den {ll\Hij\\\) = vertex. Hence a different update is 
needed, and so we implement a "small-cluster" move. 

In the small-cluster move a spin is chosen at random 
and an attempt is made to flip it. First we consider the 
case of a 1-spin. The spin in question may be connected 
to other spins (in the other integer) through bond oper- 
ators. Let us denote these spins nearest-neighbor (nn) 
spins. The nn-spins, in turn, may be connected to spins 
in the same integer as the spin we are attempting to flip. 
We call these spins next nearest neighbor (nnn) spins. 
An attempt to flip a 1-spin is always accepted, but in 
order to satisfy detailed balance with the reverse update, 
to be described below, we also need to consider all nn 
spins. If there is a nn spin connected to nnn spins, all of 
which take the value 0, then these nn spins are assigned 
values and 1 with equal probability as the original spin 
is flipped. 

The reverse move is flipping a 0-spin to a 1-spin. This 
move is accepted with probability 2~""°, where nrio is 
the number of nn spins connected to nnn spins, all of 
which necessarily assume the value . If the move is 
accepted all the nn spins are set to 0. The factor 2"""" 
corresponds to the number of ways the nn spins can be 
assigned states and 1 in the above reverse move and 
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FIG. 7: Illustration of a small-cluster update satisfying de- 
tailed balance. The upper row of spins denote one integer 
and the lower row the other integer. Spins in the 1 state are 
denoted by a filled circle, and spin spins in the 0-state be an 
empty circle. The spin marked by a square box is flipped 
between the 1-state (left side) and 0-state (right side). 



ensures that detailed balance is satisfied. An illustration 
of this move (with nno = 1) is shown in Fig. [T] The 
advantage of this update, compared to a single-spin flip, 
is that is allows 0-spins to be updated even though they 
may be connected to 1-spins. 

The two updates together ensure ergodicity, satisfy de- 
tailed balance and demonstrate that it is possible to use 
an update procedure that flips more than one spin at a 
time. However, it only works for the model Eq. ([S]), which 
does not have the prescribed integer as a ground state. 
One way to still use the small-cluster update is to allow 
only updates that do not lead to a state with an energy 
lower than the prescribed composite integer q. This we 
have implemented, but unfortunately the scaling, once 
again, is not better than simple Metropolis updates. 



E. Transverse Field 

The above algorithms are based on thermal fluctua- 
tions. However, it is also possible to modify the model 
and introduce quantum mechanical terms. Since quan- 
tum mechanical systems are able to penetrate through 
energy barriers, instead of going over energy barriers this 
method can, in certain cases be superior. We have there- 
fore added a a transverse magnetic field to the model. 



i,3 



2j 



hY^iSt + Sil (8) 



in order to compare the efficiency of the quantum and 
thermal fluctuations. First we consider a transverse field 
that is constant in time, and find the temperature and 
field strength combination that minimizes the number 
of states visited before finding the ground state. We 
implement the transverse field using a continuous time 
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FIG. 8: The average number of attempts, iV before finding a g. ^j^^ complexity, or number of MC steps, necessary to 

prime factor as a function of temperature and transverse field ^^^^^^^ ^^-^^ ^^^^^^ ^^^^ imaginary time quan- 

for a 14-spin factor. annealing method. 



algorithm^ Individual spins now fluctuate in imagi- 
nary time, S{t), and imaginary-time segments of individ- 
ual spins can be flipped. However, due to the difficulties 
explained in the above section we do not use a cluster 
update, but an imaginary-time spin segment {tq,ti} is 
flipped with probability exp(J^^ dTTAE{T). 

In Fig. [8] we show the effect of increasing the transverse 
field for the = 14 spin system. Interestingly it appears 
that, in this case, the quantum fluctuations are not more 
efficient in finding the ground state than thermal fluctu- 
ations. The smallest number of steps are found for the 
case of zero transverse field. We have also used an expo- 
nentially increasing field, hi = 1.5', which leads to more 
fluctuations in imaginary time for the higher spins, but 
we flnd that this does not improve the scaling. 



F. Quantum annealing 

While in the last subsection we considered the effect of 
a time-independent transverse field it is more common to 
gradually reduce the quantum mechanical terms during 
the simulation. A commonly used annealing method^! 
defines a Hamiltonian 

H = sH, + {s- l)H„ (9) 

where He is the classical model with the desired ground 
state. Quantum fluctuations are introduced through Hq, 
which, for Ising systems, usually is a transverse held. The 
control parameter s is slowly reduced from the initial 
value s = to the final value s = 1 as the system evolves 
according to the time-dependent Schrodinger equation 

The process is performed at a very low temperature, and, 
if the process is sufficiently slow as to be adiabatic, the 



system evolves from the ground state of Hq to the ground 
state of He, which is the solution to the problem. The 
time required to find the correct ground state with a sig- 
nificant probability is called the complexity of the prob- 
lem, and numerical simulations of small systems indicate 
that the complexity may (in some cases) increase poly- 
nomially in system sizci ^'''^^ 




FIG. 10: The 20 lowest energy levels as function of transverse 
field for an instance of the factorization problem (q=551) with 
two five bit integers. 

Real-time quantum annealing governed by the 
Schrodinger equation is limited by the size of the mini- 
mal gap between the ground state and first excited state. 
In Fig. [TUl the 20 lowest energy levels for the model de- 
fined by Eq. © are displayed. The system consist of 
two five bit integers with q — 551, and we do observe 
a minimum in the gap around h — 3. It appears that, 
at least in some random satisfiability problems, this gap 
is exponentially small for certain instances of the prob- 
lems due to a first-order phase transition j^^-^^ This effect 
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could limit the applicability of quantum annealing meth- 
ods to smaller system sizes. A recent study shows that, 
for small system sizes, real-time quantum annealing of 
the factorization model does indeed scale polynomially 
in system size.^ In order to find out whether the scaling 
persists to larger system sizes it would be necessary to 
study the scaling of the ground state energy gap with 
system sizej ^^'^^ 

Here we only use imaginary-time dynamics to study 
larger system sizes. This does not constitute the real 
time dynamics of the Schrodinger equation, but limited 
success has nevertheless been reportedi^S In order to test 
this method on the factorization problem we study the 
model 



H = 



-{s-l)hY,{St + S-). (10) 



The temperature T is set so that virtually no spin flips 
are accepted when s=l, and the strength of the trans- 
verse field is set to h = 10 T. We have determined the 
number of sweeps of the whole lattice (MC steps) neces- 
sary to find the correct ground state with a probability 
of about 30%, and display the result in Fig. |9l From 
the figure we see that the system size dependence is still 
exponential, and it appears that the standard quantum 
annealing method as applied in an imaginary-time path- 
integral simulation does not improve the scaling. 



V. DISCUSSION AND CONCLUSION 

We have developed a statistical mechanics Monte Carlo 
approach to the factorization problem. The resulting 
model is highly complex with exponentially large frus- 
trated long-range multi-spin interactions. We found that 
importance sampling is only very weakly effective in im- 
proving the convergence to the ground state. The fastest 
method we found remains exponential and beats random 
sampling by only a factor two. 

We believe that the challenge to standard statistical 
methods is threefold. First, as for other difficult opti- 
mization problems, there is a complex energy landscape 
with many low-lying states separated by high energy bar- 
riers. In this case the global minimum is required and 
therefore tempering and annealing methods that have 
been so successful in determining the low-temperature 
properties of spin glasses are of only limited use. 



Second, unlike in standard short-range spin glass mod- 
els the energy change resulting from a single spin update 
has an exponentially broad distribution, implying that 
most single-spin acceptance rates are either close to unity 
or close to zero. This gives the importance sampling the 
character of a random search. Rescaling the temperature 
and transverse fields for individual spins did not amelio- 
rate this problem. We also found that while the distri- 
bution of bond strengths was exponentially broad with 
a width increasing with system size, it was not broad 
enough to permit solution via a greedy algorithm. 

Third, as a direct consequence of the broadly dis- 
tributed energy change following a single spin fiip it fol- 
lows that there is no concept of nearness in spin space. 
Two states with nearly identical spin configurations may 
have vastly different energies. Standard importance sam- 
pling methods are based on small changes in energy, a re- 
quirement not easily satisfied for the factorization prob- 
lem. 

The statistical physics model resulting from the inte- 
ger factorization problem is a good benchmark model for 
ground state algorithms since the ground state energy is 
known by construction. We hope that this work may en- 
courage further investigations in the efficiency of possible 
cluster algorithms and annealing methods for the integer 
factorization problem. All methods considered in this 
work require 0{q'^ ) operations, to factor a composite in- 
teger q. Whether it is possible to find a stochastic method 
that scales with an exponent less than ^ remains an open 
question. Another interesting question which we have not 
addressed is the existence of a finite-temperature spin- 
glass transition. The factorization model is frustrated 
and, to some extent, disordered. These are considered 
two necessary conditions for the existence of a stable 
spin-glass phase. A freezing of the spins would certainly 
protect the ground state and further explain the difficulty 
of solving the factorization problem. 
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